% PcFingerSpelling.m
% R. Scheidt
% 5/27/04
% 

% Note that Santiago does not use wrist angles and palm arch in the current data set (Thus, sigs=19).

clear; close all

CharacterTokens='ABCDEFGHIJKLMNOPQRSTUVWXYZ';
%seq1='RJQTIKSWGLONVHXFBUZPMDAECY'; seq2='OVRZYBIMFPANSXTQHKUCWLGDJE'; seq3='TKZBOHWIRAJNMQSDXUVPLEFGYC'; seq4='CHLNXQIORDWZPAUMTBJGEYVKSF'; seq5='WTKPYOEHDMJLRAGSCUFVQNXZBI'; 
seq1='ABCDEFGHIKLMNOPQRSTUVWXY';

%RepSequence = [seq1 seq2 seq3 seq4 seq5];
RepSequence = [seq1 seq1 seq1];

%load AMA
%load Jon_FSA
load JFM_FSA_92910
CharacterSamples=[Hc];
[mm,nn]=size(CharacterSamples); NumberOfReps=mm; NumberOfSignals=19; % nn

% % remove Js and Zs
RepCount=NumberOfReps;
for rasi=RepCount:-1:1
    if(RepSequence(rasi) == 'Z')
        RepSequence(rasi) = [];
        NumberOfReps=NumberOfReps-1;
        CharacterSamples(rasi,:)=[];
    elseif(RepSequence(rasi) == 'J')
        RepSequence(rasi) = [];
        NumberOfReps=NumberOfReps-1;
        CharacterSamples(rasi,:)=[];
    end
end

Indices=zeros(26,3);
for rasi=1:26
    rep=1;
    for rasj=1:NumberOfReps
        if(RepSequence(rasj) == 'A'+(rasi-1))
            Indices(rasi,rep)=rasj;
            rep=rep+1;
        end
    end
end

% remove the mean for proper operation of the PCA
ZmAMA=zeros(NumberOfReps, NumberOfSignals);
SignalMeans=mean(CharacterSamples);
for rasi=1:NumberOfSignals
    ZmAMA(:,rasi)=CharacterSamples(:,rasi)- SignalMeans(rasi);
end

%  work on the means - note that these are already zero-mean!
NumberOfLetters = 26;
Letters = zeros(NumberOfLetters, NumberOfSignals);
Jindex=10;
Zindex=26;
for rasi=1:NumberOfLetters
    if ((rasi==Jindex)|(rasi==Zindex)) % if rasi points to a J or an Z
        continue;
    end
    Letters(rasi,:) = mean(ZmAMA(Indices(rasi,:),:));
end
Letters(Zindex,:)=[];
Letters(Jindex,:)=[];

Sigma=Letters'*Letters;  % this should be a signal x signal matrix

[V,D] = eig(Sigma);
[D,I]=sort(diag(D));% find the column number with the greatest eigenvalue in D... (I think it is the last element)
D=flipud(D);  % equivalent to diag(S) from svd(Sigma)
I=flipud(I);
W=V(:,I);

% select only 1st two (or five, or N_PC) eigenvectors (principal components)
NumberOfPCs=2;
PC=zeros(size(W));
PC(:,1:NumberOfPCs)=W(:,1:NumberOfPCs);

% obtain the expansion coefficients for each manual letter (20 signal vector) 
%    Rotate along primary PCs, throwing away higher order info
coeff=Letters(:,1:NumberOfSignals)*PC; % PC coefficients for each manual character - note: (R'*X')' is equivalent to X*R
ProjLetters = coeff*PC';% RE-rotate reduced dimensional representation into the full-dimensional space
Residuals=Letters-ProjLetters;% Now create error vector (residuals in the original full Dimensional space)

TempVariance=zeros(1,NumberOfSignals);
for rasj=1:NumberOfSignals
    TempVariance(rasj) = Letters(:,rasj)'*Letters(:,rasj); % implicit summation over the entire set of (original) characters 
end
TotalVariance=sum(TempVariance);  % sum over the signals
TempResidVariance=zeros(1,NumberOfSignals);
for rasj=1:NumberOfSignals
    TempResidVariance(rasj) = Residuals(:,rasj)'*Residuals(:,rasj); % implicit summation over the entire set of (rotated) characters 
end
ResidualVariance=sum(TempResidVariance);
VAFOverSet = (1-ResidualVariance./TotalVariance)*100;

Results.NumberOfPCs = NumberOfPCs;
Results.TotalMeanLetterVariance = TotalVariance;
Results.ResidualMeanLetterVariance = ResidualVariance;
Results.OverallMeanLetterVAF = VAFOverSet;

% plot average Letters
figure; hold on
AMAcharacter = 'A'-1;
LetterIndex = 0;
for rasi=1:NumberOfLetters
    AMAcharacter=AMAcharacter+1;
    if ((rasi==Jindex)|(rasi==Zindex)) % if rasi points to a J or an Z
        continue;
    end
    LetterIndex=LetterIndex+1;
    plot(coeff(LetterIndex,1),coeff(LetterIndex,2),'w.');
    text(coeff(LetterIndex,1),coeff(LetterIndex,2),char(AMAcharacter));
end

% plot all reps
RepCoeff=ZmAMA*PC; % PC coefficients for each manual character - note: (R'*X')' is equivalent to X*R
figure; hold on
for rasi=1:NumberOfLetters
    if ((rasi==Jindex)|(rasi==Zindex)) % if rasi points to a J or an I
        continue;
    end
    plot(RepCoeff(Indices(rasi,:),1),RepCoeff(Indices(rasi,:),2),'w.');
    text(RepCoeff(Indices(rasi,:),1),RepCoeff(Indices(rasi,:),2),RepSequence(Indices(rasi,:))');
end


% -------------------------------------------------------------------------------------------------------------
% Bisection Method (bi-by-eye)
% -------------------------------------------------------------------------------------------------------------
% CharacterTokens='ABCDEFGHIJKLMNOPQRSTUVWXYZ';
% set1=[1 2 3 5 8 9 12 13 14 15 20 22];
% set2=[4 6 7 11 16 17 18 19 21 23 24 25];
% set1=[3 4 8 9 12 22];
% set2=[1 2 5 6 7 11 13 14 15 16 17 18 19 20 21 23 24 25];
% set1=[1 2 3 5 13 14 15 19 20];  % A B C E M N O S T
% set2=[4 6 7 8 9 11 12 16 17 18 21 22 23 24 25];
set1=[2 3 6 22 23];  % B C F V W
set2=[1 4 5 7 8 9 11 12 13 14 15 16 17 18 19 20 21 24 25];
set3=[10 26];  % throw J and Z away

set1=[2 3 6 22 23];  % B C F V W     VAFOverSet1 = 79.2653
set2=[1 8 13 19 24]; % A H M S X     Set2VAFOverSet1 = 26.2234

% set1=[1 2 3 22 23];  % A B C F V W   VAFOverSet1 = 80.9173
% set1=[1 2 3 (22-1) (23-1)];  % A B C V W   VAFOverSet1 = 81.5390 correct for missing 10 and 26
set1=[1 2 3 (22-1) (23-1)];  % A B C V W   VAFOverSet1 = 81.5390 correct for missing 10 and 26
% set1=[1 2 (22-1) (23-1) 25-1];  % A B V W Y   VAFOverSet1 = 76.4845
% set2=[5 7 14 17 20]; % E G N Q T     Set2VAFOverSet1 = 16.1665
% set2=[5 7 14-1 17-1 20-1]; % E G N Q T     Set2VAFOverSet1 = 21.4726
% set2=[5 8 11 20 24]; % E H K T X     Set2VAFOverSet1 = 15.7725
% set2=[5 8 (11-1) (20-1) (24-1)]; % E H K T X     Set2VAFOverSet1 = 22.5109
set2=[5 8 (12-1) (21-1) (25-1)]; % E H L U Y     Set2VAFOverSet1 = 15.7725
% set3=[4 12 13 16 18]; % D L M P R    Set3VAFOverSet1 = 19.5900
% set3=[15-1 19-1 20-1 21-1 25-1]; % O S T U Y  Set3VAFOverSet1 = 25.4547
% set3=[6 9 12-1 13-1 16-1]; % F I L M P   Set3VAFOverSet1 = 30.6562
% set3=[6 9 14-1 19-1 20-1]; % F I N S T  Set3VAFOverSet1 = 39.0229
set3=[6 14-1 19-1 20-1 21-1]; % F I N T U Set3VAFOverSet1 = 40.1808
% set3=[4 12-1 15-1 19-1 24-1]; % D L O S X  Set3VAFOverSet1 = 29.8932

set1=[1 2 3 22-1 23-1];  % A B C V W       VAFOverSet1 = 81.5390 correct for missing 10 and 26
set2=[6 14-1 19-1 20-1 21-1]; % F I N T U  Set2VAFOverSet1 = 40.1808
set3=[4 5 13-1 24-1 25-1]; % D E M X Y     Set3VAFOverSet1 = 19.8509
set1char = 'ABCVW'; set2char = 'FINTU'; set3char = 'DEMXY';


if (1)   % set 1
    SelectedTrials = reshape(Indices(set1,:),1,[]);
%     SelectedLetters = ZmAMA(SelectedTrials,:);
    SelectedLetters = Letters(set1,:);
    Sigma=SelectedLetters'*SelectedLetters;
    
    [V,D] = eig(Sigma); [D,I]=sort(diag(D)); D=flipud(D); I=flipud(I); W=V(:,I);
    NumberOfSet1PCs = 4; PC=zeros(size(W)); PC(:,1:NumberOfSet1PCs)=W(:,1:NumberOfSet1PCs);
    
    RAS_PC2 = PC;  % Xiaolin add code here!
    
    PC_From_RAS = RAS_PC2(:,1:4);
    save JFM_92910_PC2 PC_From_RAS
    
    
    coeff=SelectedLetters*PC; % PC coefficients for each manual character - note: (R'*X')' is equivalent to X*R
    
    ProjLetters = coeff*PC';% RE-rotate reduced dimensional representation into the full-dimensional space
    Residuals=SelectedLetters-ProjLetters;% Now create error vector (residuals in the original full Dimensional space)
    
    TempVariance=zeros(1,NumberOfSignals);
    for rasj=1:NumberOfSignals
        TempVariance(rasj) = SelectedLetters(:,rasj)'*SelectedLetters(:,rasj); % implicit summation over the (original) characters 
    end
    TotalVariance=sum(TempVariance);  % sum over the signals
    TempResidVariance=zeros(1,NumberOfSignals);
    for rasj=1:NumberOfSignals
        TempResidVariance(rasj) = Residuals(:,rasj)'*Residuals(:,rasj); % implicit summation over the (original) characters 
    end
    ResidualVariance=sum(TempResidVariance);
    
    % The following were added by Xiaolin
    ResidualVariance
    TotalVariance
    % Added by Xiaolin - End
      
    VAFOverSet1 = (1-ResidualVariance./TotalVariance)*100
    Results.NumberOfSet1PCs = NumberOfSet1PCs;
    Results.Set1PCs = PC;
    Results.TotalSet1LetterVariance = TotalVariance;
    Results.ResidualSet1LetterVariance = ResidualVariance;
    Results.OverallSet1LetterVAF = VAFOverSet1;
 
    Set2SelectedLetters = Letters(set2,:);
    Set2coeff=Set2SelectedLetters*PC; % PC coefficients for each manual character - note: (R'*X')' is equivalent to X*R
    ProjSet2Letters = Set2coeff*PC';% RE-rotate reduced dimensional representation into the full-dimensional space
    Set2Residuals=Set2SelectedLetters-ProjSet2Letters;% Now create error vector (residuals in the original full Dimensional space)
    TempVariance=zeros(1,NumberOfSignals);
    for rasj=1:NumberOfSignals
        TempVariance(rasj) = Set2SelectedLetters(:,rasj)'*Set2SelectedLetters(:,rasj); % implicit summation over the (original) characters 
    end
    Set2TotalVariance=sum(TempVariance);  % sum over the signals
    TempResidVariance=zeros(1,NumberOfSignals);
    for rasj=1:NumberOfSignals
        TempResidVariance(rasj) = Set2Residuals(:,rasj)'*Set2Residuals(:,rasj); % implicit summation over the (original) characters 
    end
    Set2ResidualVariance=sum(TempResidVariance);
    Set2VAFOverSet1 = (1-Set2ResidualVariance./Set2TotalVariance)*100

    Set3SelectedLetters = Letters(set3,:);
    Set3coeff=Set3SelectedLetters*PC; % PC coefficients for each manual character - note: (R'*X')' is equivalent to X*R
    
    figure; hold on; title('set 1 PC space')% plot average Letters
    LetterIndex = 0;
    for rasi=1:length(set1)
        if ((rasi==Jindex)|(rasi==Zindex)) % if rasi points to a J or an Z
            continue;
        end
        LetterIndex=LetterIndex+1;
        plot(coeff(LetterIndex,1),coeff(LetterIndex,2),'w.');
        AMAcharacter=set1char(rasi);
        text(coeff(LetterIndex,1),coeff(LetterIndex,2),char(AMAcharacter),'FontSize',20);
        AMAcharacter=set2char(rasi);
        text(Set2coeff(LetterIndex,1),Set2coeff(LetterIndex,2),char(AMAcharacter),'FontSize',12);
        AMAcharacter=set3char(rasi);
        text(Set3coeff(LetterIndex,1),Set3coeff(LetterIndex,2),char(AMAcharacter),'FontSize',12);
    end

    Set2Sigma=Set2SelectedLetters'*Set2SelectedLetters;
    [V,D] = eig(Set2Sigma); [D,I]=sort(diag(D)); D=flipud(D); I=flipud(I); W=V(:,I);
    NumberOfSet1PCs = 2; PC=zeros(size(W)); PC(:,1:NumberOfSet1PCs)=W(:,1:NumberOfSet1PCs);
    Set2coeff=Set2SelectedLetters*PC; % PC coefficients for each manual character - note: (R'*X')' is equivalent to X*R
    ProjSet2Letters = Set2coeff*PC';% RE-rotate reduced dimensional representation into the full-dimensional space
    Set2Residuals=Set2SelectedLetters-ProjSet2Letters;% Now create error vector (residuals in the original full Dimensional space)
    TempVariance=zeros(1,NumberOfSignals);
    for rasj=1:NumberOfSignals
        TempVariance(rasj) = Set2SelectedLetters(:,rasj)'*Set2SelectedLetters(:,rasj); % implicit summation over the (original) characters 
    end
    Set2TotalVariance=sum(TempVariance);  % sum over the signals
    TempResidVariance=zeros(1,NumberOfSignals);
    for rasj=1:NumberOfSignals
        TempResidVariance(rasj) = Set2Residuals(:,rasj)'*Set2Residuals(:,rasj); % implicit summation over the (original) characters 
    end
    Set2ResidualVariance=sum(TempResidVariance);
    Set2VAFOverSet2 = (1-Set2ResidualVariance./Set2TotalVariance)*100

    SelectedLetters = Letters(set1,:);
    coeff=SelectedLetters*PC; % PC coefficients for each manual character - note: (R'*X')' is equivalent to X*R
    Set3SelectedLetters = Letters(set3,:);
    Set3coeff=Set3SelectedLetters*PC; % PC coefficients for each manual character - note: (R'*X')' is equivalent to X*R

    figure; hold on; title('set 2 PC space')% plot average Letters
    LetterIndex = 0;
    for rasi=1:length(set2)
        AMAcharacter=set2char(rasi);
        if ((rasi==Jindex)|(rasi==Zindex)) % if rasi points to a J or an I
            continue;
        end
        LetterIndex=LetterIndex+1;
        plot(Set2coeff(LetterIndex,1),Set2coeff(LetterIndex,2),'w.');
        AMAcharacter=set1char(rasi);
        text(coeff(LetterIndex,1),coeff(LetterIndex,2),char(AMAcharacter),'FontSize',12);
        AMAcharacter=set2char(rasi);
        text(Set2coeff(LetterIndex,1),Set2coeff(LetterIndex,2),char(AMAcharacter),'FontSize',20);
        AMAcharacter=set3char(rasi);
        text(Set3coeff(LetterIndex,1),Set3coeff(LetterIndex,2),char(AMAcharacter),'FontSize',12);
    end

    
    Set3SelectedLetters = Letters(set3,:);
    Set3coeff=Set3SelectedLetters*PC; % PC coefficients for each manual character - note: (R'*X')' is equivalent to X*R
    ProjSet3Letters = Set3coeff*PC';% RE-rotate reduced dimensional representation into the full-dimensional space
    Set3Residuals=Set3SelectedLetters-ProjSet3Letters;% Now create error vector (residuals in the original full Dimensional space)
    TempVariance=zeros(1,NumberOfSignals);
    for rasj=1:NumberOfSignals
        TempVariance(rasj) = Set3SelectedLetters(:,rasj)'*Set3SelectedLetters(:,rasj); % implicit summation over the (original) characters 
    end
    Set3TotalVariance=sum(TempVariance);  % sum over the signals
    TempResidVariance=zeros(1,NumberOfSignals);
    for rasj=1:NumberOfSignals
         TempResidVariance(rasj) = Set3Residuals(:,rasj)'*Set3Residuals(:,rasj); % implicit summation over the (original) characters 
    end
    Set3ResidualVariance=sum(TempResidVariance);
    Set3VAFOverSet1 = (1-Set3ResidualVariance./Set3TotalVariance)*100
    
    
    
    
    
    if(1)  % do this over reps
        RepCoeff=ZmAMA*PC; % project all characters onto this plane
        ProjZmAMA = RepCoeff*PC'; % then RE-rotate reduced dimensional representation into the full-dimensional space
        TotalCharVariance=zeros(NumberOfLetters,1);
        PCCharVariance=zeros(NumberOfLetters,1);
        LetterIndex = 0;
        for rasi=1:NumberOfLetters
            LetterIndex=LetterIndex+1;
            if ((rasi==Jindex)|(rasi==Zindex)) % if rasi points to a J or an I
                continue;
            end
            
            TempCharVariance=zeros(1,NumberOfSignals);
            for rasj=1:NumberOfSignals
                CharacterSignal = ZmAMA(Indices(LetterIndex,:),rasj);
                CharacterSignal= CharacterSignal-mean(CharacterSignal);  % subtract the mean of the signal of interest over the repetitions of the character
                TempCharVariance(rasj) = CharacterSignal'*CharacterSignal; % implicit summation over the 5 repetitions
            end
            TotalCharVariance(rasi)=sum(TempCharVariance);
            
            TempPCVariance=zeros(1,NumberOfPCs);
            for rasj=1:NumberOfSignals
                CharacterSignal = ProjZmAMA(Indices(LetterIndex,:),rasj);
                CharacterSignal= CharacterSignal-mean(CharacterSignal);  % subtract the mean of the signal of interest over the repetitions of the character
                TempPCVariance(rasj) = CharacterSignal'*CharacterSignal; % implicit summation over the 5 repetitions
            end
            PCCharVariance(rasi)=sum(TempPCVariance);
        end
%         disp(CharacterTokens(LetterIndex));
        VAFByCharacter = (PCCharVariance./TotalCharVariance)*100;
        Results.Subset1.list = set1;
        Results.Subset1.VAFByCharacter = VAFByCharacter;
    end
end

% end

Results

break


if (1)   % set 2
    SelectedTrials = reshape(Indices(set2,:),1,[]);
    SelectedLetters = ZmAMA(SelectedTrials,:);
    Sigma=SelectedLetters'*SelectedLetters;
    
    [V,D] = eig(Sigma); [D,I]=sort(diag(D)); D=flipud(D); I=flipud(I); W=V(:,I);
    NumberOfSet2PCs = 2; PC=zeros(size(W)); PC(:,1:NumberOfSet2PCs)=W(:,1:NumberOfSet2PCs);
    
    coeff=SelectedLetters*PC; % PC coefficients for each manual character - note: (R'*X')' is equivalent to X*R
    ProjLetters = coeff*PC';% RE-rotate reduced dimensional representation into the full-dimensional space
    Residuals=SelectedLetters-ProjLetters;% Now create error vector (residuals in the original full Dimensional space)
    
    TempVariance=zeros(1,NumberOfSignals);
    for rasj=1:NumberOfSignals
        TempVariance(rasj) = SelectedLetters(:,rasj)'*SelectedLetters(:,rasj); % implicit summation over the (original) characters 
    end
    TotalVariance=sum(TempVariance);  % sum over the signals
    TempResidVariance=zeros(1,NumberOfSignals);
    for rasj=1:NumberOfSignals
        TempResidVariance(rasj) = Residuals(:,rasj)'*Residuals(:,rasj); % implicit summation over the (original) characters 
    end
    ResidualVariance=sum(TempResidVariance);
    VAFOverSet2 = (1-ResidualVariance./TotalVariance)*100;
    
    Results.NumberOfSet2PCs = NumberOfSet2PCs;
    Results.Set2PCs = PC;
    Results.TotalSet2LetterVariance = TotalVariance;
    Results.ResidualSet2LetterVariance = ResidualVariance;
    Results.OverallSet2LetterVAF = VAFOverSet2;
    
    if(1)  % do this over reps
        RepCoeff=ZmAMA*PC; % project all characters onto this plane
        ProjZmAMA = RepCoeff*PC'; % then RE-rotate reduced dimensional representation into the full-dimensional space
        TotalCharVariance=zeros(NumberOfLetters,1);
        PCCharVariance=zeros(NumberOfLetters,1);
        LetterIndex = 0;
        for rasi=1:NumberOfLetters
            LetterIndex=LetterIndex+1;
            if ((rasi==Jindex)|(rasi==Zindex)) % if rasi points to a J or an I
                continue;
            end
            
            TempCharVariance=zeros(1,NumberOfSignals);
            for rasj=1:NumberOfSignals
                CharacterSignal = ZmAMA(Indices(LetterIndex,:),rasj);
                CharacterSignal= CharacterSignal-mean(CharacterSignal);  % subtract the mean of the signal of interest over the repetitions of the character
                TempCharVariance(rasj) = CharacterSignal'*CharacterSignal; % implicit summation over the 5 repetitions
            end
            TotalCharVariance(rasi)=sum(TempCharVariance);
            
            TempPCVariance=zeros(1,NumberOfPCs);
            for rasj=1:NumberOfSignals
                CharacterSignal = ProjZmAMA(Indices(LetterIndex,:),rasj);
                CharacterSignal= CharacterSignal-mean(CharacterSignal);  % subtract the mean of the signal of interest over the repetitions of the character
                TempPCVariance(rasj) = CharacterSignal'*CharacterSignal; % implicit summation over the 5 repetitions
            end
            PCCharVariance(rasi)=sum(TempPCVariance);
        end
%         disp(CharacterTokens(LetterIndex));
        VAFByCharacter = (PCCharVariance./TotalCharVariance)*100;
        Results.Subset2.list = set2;
        Results.Subset2.VAFByCharacter = VAFByCharacter;
    end
    
end

Results
